import numpy as np
import pandas as pd
import matlab.engine

eng = matlab.engine.start_matlab()
import scipy.io as scio


def cal_pz(H, errolist, k):
    for i in range(k):
        Hs = H[:, i]


def test_pz(Hs, rate,flag):
    if flag ==1:
        print(flag)
    # 计算单个直方图分别在0-254喝255-512中大于0的点的个数
    # Left = np.sum(Hs[0:254] > 0)  # -255-0
    # Right = np.sum(Hs[255:512] > 0)  # 256-512
    # pp = np.zeros((1, Left))#记录1-256中非0的元素
    # nps = np.zeros((1, Right))#记录-255-0中非0元素的坐标
    # pz = np.zeros((1, 255 - len(pp)))#记录1-256中0的元素
    # nz = np.zeros((1, 255 - len(nps)))#记录-255-0中0元素的坐标
    nps = np.where(Hs[0:255] != 0)[0]  # 记录1-256中非0的元素
    nz = np.where(Hs[0:255] == 0)[0]  # 记录1-256中0的元素
    pp = np.where(Hs[255:512] != 0)[0]  # 记录-255-0中非0元素的坐标
    pz = np.where(Hs[255:512] == 0)[0]  # 记录-255-0中0元素的坐标
    nps = -np.sort(-nps)
    nz = -np.sort(-nz)
    print('nps',len(nps))
    print('pp',len(pp))
    maxT = 0  # 最近的0值点
    minT = 0
    if 255 - nz[0] < pz[0]:
        maxT = 254 - nz[0]
    else:
        maxT = pz[0]
    s = Hs[255]
    # 找到第一个满足分配容量的点
    for i in range(maxT):
        if s >= rate:
            break
        if i > 0:
            s = s + Hs[255 + i] + Hs[255 - i]
            minT = i
    if minT == maxT + 1:
        P = 0
        Z = 0
        return P, Z
    tmp = 0

    while Hs[tmp] == 0 :
        if tmp<510:
            tmp += 1
        else:
            break
    print('asdas',tmp)
    nz = np.arange(tmp, 0, -1)

    tmp = 510
    while Hs[tmp] == 0:
        if tmp>0:
            tmp -= 1
        else:
            break
    pz = np.arange(tmp + 1, 511, 1)
    tmd = 100000000

    P = []
    Z = []
    Rate = 0
    D = 0
    for t in range(minT, maxT):
        R = 0
        D = 0
        i = t
        tP = []
        tZ = []
        while R < rate:
            R = R + Hs[i + 255]
            if i >= 0:
                tP.append(pp[i] + 255)
                if t - i > len(pz):
                    P = [254, -254]
                    Z = [255, -255]
                    return
                tZ.append(pz[t - i])
                i = -i
            else:
                tP.append(nps[-i - 1])
                if t + i > len(nz):
                    P = [254, -254]
                    Z = [255, -255]
                    return
                # print(t + i+1)
                # print(nz[t + i],nz[t + i+1],nz[t + i-1])
                tZ.append(nz[t + i+1])
                i = -i - 1

        R, D = rateAndDist(tP, tZ, Hs)
        if D < tmd:
            tmd = D
            P = tP
            Z = tZ
            Rate = R
    print('p,z,r,d', np.array(P) - 255, np.array(Z) - 255, Rate, tmd)
    return P, Z


def rateAndDist(Ps, Zs, H):
    # print('p', Ps)
    m = len(Ps)
    n = len(H)
    Ps = np.array(Ps, dtype=int)
    Zs = np.array(Zs, dtype=int)
    C = np.zeros((1, n), dtype='double')
    PK = np.zeros((1, n), dtype='double')
    sig = np.sign(Zs - Ps)
    PK[0][Ps] = sig
    # print(PK,Ps,sig,Zs)
    for i in range(m):
        # print(Ps[i] + sig[i], Zs[i] - sig[i] + 1, sig[i])
        tem = np.arange(Ps[i] + sig[i], Zs[i] - sig[i] + 1, sig[i])
        C[0, tem] += sig[i]
    rates = sum(H[Ps])
    # print(C)
    # print(rates)
    H = H.reshape(511, -1)
    dists = 0.5 * (np.sum(np.dot((np.square(C) + np.square(C + PK)), H)))
    # print(rates, dists)
    return rates, dists


def initation(numV, maxL, H, PX, PZ):
    '''

    :param numV: 种群数量
    :param maxL: 最大维度
    :param H: 直方图
    :param PX: 峰值点
    :param PZ: 0值点
    :return:
    '''
    V = np.zeros((numV, maxL * 3))
    tpos = np.argsort(-H)  # 降序排序索引
    th = H[tpos]  # 降序排序直方图
    numz = np.sum(th == 0)  # 0值点个数
    nump = np.sum(th > 18)  # 峰值点维度个数
    if nump < maxL:
        nump = maxL
    # 初始化种子
    pSet = tpos[0:nump]
    zSet = tpos[len(H) - numz:len(H)]
    tempP = pSet.copy()  # 初始化Ni的经验种子
    tempZ = zSet.copy()
    V[0, maxL:2 * maxL] = pSet[0:maxL]
    for i in range(maxL):
        # 找到与峰值点最近的一个0值点,将其塞入种群当中
        index = np.argmin(tempZ - V[0, maxL + i])
        V[0, 2 * maxL + i] = tempZ[index]
        tempZ[index] = 1024
    # 正叙从000000-111111
    for i in range(maxL):
        V[i, 0:i + 1] = np.ones((1, i + 1))
        V[i, maxL: maxL * 3] = V[0, maxL: maxL * 3]
    # 找到pSet中不属于Px的元素，并正序取到每一个值
    tempP = np.setdiff1d(pSet, PX)
    tempZ = np.setdiff1d(zSet, PZ)
    npx = len(PX)
    # 从maxL-2*maxL,
    for i in range(maxL, 2 * maxL):
        V[i, 0:npx] = np.ones((1, npx))
        V[i, npx:maxL] = (np.random.rand(maxL - npx) > 0.5).astype(int)
        V[i, maxL:maxL + npx] = PX
        V[i, 2 * maxL:2 * maxL + npx] = PZ
        # 从剩下的maxL-npx随机从tempP,temzZ中取填补上
        a = np.arange(len(tempP))
        np.random.shuffle(a)
        b = np.arange(len(tempZ))
        np.random.shuffle(b)
        for j in range(maxL - npx):
            V[i, maxL + npx + j] = tempP[a[j]]
            V[i, 2 * maxL + npx + j] = tempZ[b[j]]
    # 剩下的2*maxl-numv将pset,vset随机插入
    for i in range(2 * maxL, numV):
        a = np.arange(len(pSet))
        np.random.shuffle(a)
        for j in range(maxL):
            V[i, maxL + j] = pSet[a[j]]
        a = np.arange(len(zSet))
        np.random.shuffle(a)
        for j in range(maxL, 2 * maxL):
            V[i, maxL + j] = zSet[a[j]]
        V[i, 0:maxL] = (np.random.rand(maxL) > 0.1).astype(int)
    return V, pSet, zSet


def fitness(numV_s, Vs, maxL, minL, H):
    '''
    :param numV_s: 种群数量
    :param Vs: 初始化种群编码
    :param maxL: 最大维度
    :param minL: 最少容量
    :param H: 直方图
    :return:
    '''
    # print(numV_s)
    fit = np.zeros(numV_s)
    valid = np.ones(numV_s)
    # 将所有取1的点计算其失真加入到fit中
    for i in range(numV_s):
        thum = int(np.sum(Vs[i, 0:maxL]))
        Px = np.ones(int(thum))
        Zx = Px.copy()
        for j in range(maxL - 1, -1, -1):

            if Vs[i, j]:
                Px[thum - 1] = Vs[i, maxL + j]

                Zx[thum - 1] = Vs[i, 2 * maxL + j]
                thum -= 1

        [R, D] = rateAndDist(Px, Zx, H)
        fit[i] = D
        if R < minL + 18 * thum:
            valid[i] = 0
    maxnum = np.max(fit)
    for i in range(numV_s):
        if valid[i]:
            fit[i] = maxnum - fit[i] + 1
        else:
            fit[i] = 1
    # 归一化
    Afit = fit / np.sum(fit)
    Afit=np.cumsum(Afit)

    return Afit, fit


def selection(Afit):
    '''
    :param Afit: 归一化概率
    :return: 返回两个随机适度体
    '''
    c1 = selectone(Afit, 0, len(Afit))
    c2 = c1
    while c1 == c2:
        c2 = selectone(Afit, 0, len(Afit))
    return c1, c2


def selectone(AFit, lb, ub):
    '''
    :param AFit:归一化适度概率
    :param lb:上限
    :param ub:下限
    :param m:随机概率
    :return: 返回随机适度个体
    '''
    # idx =0
    m =np.random.rand(1)[0]

    if lb == ub:
        idx = lb
    else:
        mid = np.floor((lb + ub-1)/2)
        mid = int(mid)
        if m < AFit[mid]:
            idx = selectone(AFit, lb, mid)
        else:
            print(9999999)
            idx = selectone(AFit, mid+1, ub)
    return idx


def crossover(fit1, fit2, pSet, zSet, maxL):
    k = np.floor(np.random.rand(1)[0] * 2 * maxL)
    k = int(k)
    c = fit1.copy()
    c[maxL + k: 3 * maxL] = fit2[maxL + k: 3 * maxL]
    fit2[maxL + k: 3 * maxL] = fit1[maxL + k: 3 * maxL]
    fit1 = c.copy()

    return fit1, fit2


def mutation(fit1, fit2, pSet, zSet, maxL):
    # k = np.floor(np.random.rand(1)[0] * 3 * maxL)
    i = np.random.rand(maxL)[0]
    i =int(i)
    if fit1[i] == 0:
        fit1[i] = 1
    else:
        if fit1[i] == 1:
            fit1[i] = 0
    if fit2[i] == 0:
        fit2[i] = 1
    else:
        if fit2[i] == 1:
            fit2[i] = 0
    return fit1, fit2


def Ga_method(H, Gs, pcs, pms, numVs, maxL, minR):
    g = 0  # 代数
    Vs = []  # 种群
    pSet = []  # 峰值点集
    zSet = []  # 零点集
    Ps = []
    Zs = []
    bestDist = np.zeros( Gs)  # 每一代的最小失真
    fit = []  # 个体适应值
    best1 = np.zeros((1, maxL * 3))  # 每一代保留的最优个体
    best2 = best1.copy()  # 每一代保留的次优个体
    PX, ZX = test_pz(H, minR)
    Vs, pSet, zSet = initation(numVs, maxL, H, PX, ZX)
    nv = Vs.copy()
    # AFit, fit = fitness(numVs, Vs, maxL, minR, H)

    while g < Gs-1:
        print('第{}轮'.format(g))
        AFit, fit = fitness(numVs, Vs, maxL, minR, H)
        tempFit = fit.copy()  # 保存两个最优个体
        ind = np.argmax(tempFit)
        tempFit[ind] = 1
        best1 = Vs[ind, :].copy()
        ind = np.argmax(tempFit);
        best2 = Vs[ind, :].copy()
        j = 0
        while j < numVs - 2:
            aidx, bidx = selection(AFit)  # 选择两个个体
            a = Vs[aidx, :]
            b = Vs[bidx, :]
            if np.random.rand(1)[0] < pcs:  # 交叉
                a, b = crossover(a, b, pSet, zSet, maxL)
            if np.random.rand(1)[0] < pms:  # 变异
                a, b = mutation(a, b, pSet, zSet, maxL)
            nv[j, :] = a
            nv[j + 1, :] = b
            j = j + 2
        nv[numVs - 2, :] = best1  # 保留最优两个个体到子代
        nv[numVs-1, :] = best2
        Vs = nv
        tnum = int(np.sum(best1[1:maxL]))  # 记录上一代最优个体的失真
        Ps = np.ones(int(tnum))
        Zs = Ps.copy()
        for i in range(maxL - 1, 0, -1):
            if best1[i] and tnum > 0:
                print('ss',tnum)
                # print('bbb',best1(maxL + i))
                Ps[tnum - 1] = best1[maxL + i]
                Zs[tnum - 1] = best1[maxL * 2 + i]
                tnum = tnum - 1
        rates, dists = rateAndDist(Ps, Zs, H)
        g = g + 1
        bestDist[g] = dists
    # print('caonima',numVs)
    AFit, fit = fitness(numVs, Vs, maxL, minR, H) # 输出最优个体
    k = np.argmax(fit)
    tnum = int(np.sum(Vs[k, 1:maxL]))
    Ps = np.ones(int(tnum))
    Zs = Ps.copy()
    for i in range(maxL - 1, 0, -1):
        if Vs[k, i] and tnum > 0:
            Ps[tnum - 1] = Vs[k, maxL + i]
            Zs[tnum - 1] = Vs[k, maxL * 2 + i]
            tnum = tnum - 1
    rates, dists = rateAndDist(Ps, Zs, H)
    return Ps, Zs, rates, dists, bestDist


if __name__ == '__main__':
    # data = pd.read_csv('./H.csv').values
    # data = data[:,1]
    # print(data[255])
    err = eng.readdata('sorted_errors1.mat')
    # print(err)
    err = np.array(err['sorted_errors'][0])
    data = eng.readdata('H2.mat')
    data = np.array(data['H'][0])

    # print(data.shape)
    P,Z=test_pz(data,6554)
    # P = np.array([3, -3, 2, -2, 1, -1, 0]) + 255
    # Z = np.array([15, -15, 16, -16, 17, -17, 18]) + 255
    # P = np.array([3, ]) + 255
    # Z = np.array([15]) + 255
    # print(rateAndDist(P,Z,data))
    V = initation(100, 14, data, P, Z)
    file_names = 'Vs.mat'
    scio.savemat(file_names, {'Vs': V})
    maxl = 25
    G = 100
    pc = 0.8
    pm = 0.2
    numV = 100
    if maxl > np.sum(data > 18):
        maxl = np.sum(data > 18)
    Ps, Zs, rates, dists, bestDist=Ga_method(data, G, pc, pm, numV, maxl, 6554)
    print(Ps, Zs, rates, dists, bestDist)
